Majorana fermions on the Abrikosov flux lattice in a p x + ip y superconductor and 

thermal conductivity in superclean regime 
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We show that periodic lattice ol Abrikosov vortices in chiral p wave superconductor in general 
supports fermionic states with zero energy provided the intervortex distance is smaller than the 
critical one. The zero modes appear at the intersection with the Fermi level ol electronic magnetic 
Bloch bands formed by the overlapping vortex core states. The Bloch bands are robust against 
lattice disorder induced by fluctuations ol vortex positions and can transmit the energy flow across 
the lattice. The hallmark ol zero modes on Bloch bands in electronic heat conductivity is discussed. 

PACS numbers: 



Topological 2+1 dimensional Fermi systems are one 
of the most intriguing topics in the field of condensed 
matter physics. An intense investigation of such systems 
has started from the theory of integer quantum Hall ef- 
fect when it was shown that the transverse conductivity 
is proportional to the discrete valued topological invari- 
ant of the ground stated related to the first Chern num- 
ber of the Berry phase gauge field in the Brillpuin zone. 
An analogous topologically nontrivial ground state was 
found in superfluid A phase of 3 He films^. In this system 
the non-trivial chiral p x + ip y structure of superfuid or- 
der parameter corresponding to the Cooper pairing with 
angular momentum L z = ±1 allows for the existence of 
quantum Hall effect in the absence of magnetic field. The 
same state is suggested to realize in layered triplet p wave 
superconductor Sr 2 RuO^. 

Many of the exotic properties of chiral p x + ip y super- 
conductors and Fermi superfluids are determined by the 
interplay of the ground state topology and the proper- 
ties of fermionic bound states modified in the vicinity of 
topological defects in order parameter distribution. In 
particular the fermionic sectors of 3 He A and Sr2Ru04 
contain zero energy states localized near domain walls 
and solitonsi, boundaries^ and quantized vortices^. The 
zero energy fermionic modes can be described in terms 
of the self-conjugated Majorana fermions which were 
theoretically predicted to appear in several other two- 
dimensional systems such as the fractional quantum Hall 
liquid at filling 5/2&, heterostructures of topological in- 
sulators and superconductors^, and possibly certain Iri- 
dates which effectively realize the Kitaev honeycomb 
models. 

An appealing possibility offered by the nontrivial 
structure of fermionic spectrum in the vortex phase of 
chiral p x + ip y superconductors is the realization of 
quantum matter with exotic non-Abelian quasiparticle 
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In this case the non-Abelian anyons are 
presented by vortex excitations supporting zero-energy 
Majorana fermions residing inside their cores. That is 
the spectrum of vortex core fermions is given by 

e = uj(n + "f) (1) 

where n is integer number, 7=1/2 for s wave£ and 7 = 



(a)_ 

! Quasiparticle unit cell 



(b) 






Ka/n 



FIG. 1: (a) The vortex lattice unit cell and corresponding 
quasiparticle unite cell containing the nodes marked by 1 and 
2 in square vortex lattice. The vortex positions are marked 
with red filled circles which form the 2D Brave lattice with 
the basis (a, b). (b) The 3D plot ol magnetic Bloch band 
e +(K) given by J52]) x = and 7 a = 7 a = 0.14. (c) The 
Bloch bands e ±{K x ,K y = 0) for 7 a = I b = 1.4; 2.8; V2tt 
shown by red dashed, blue dash-dotted and black solid lines 
correspondingly, (d) The contour plot ol eo+(K) for \ — 7r/4 
and I a = lb = 4.3. The plots (b,c,d) correspond for the 
square vortex lattice. 



for p x + ip y wave^ superconductors. Thus in topolog- 
ically non-trivial superconductors the spectrum of vor- 
tex core fermions ([1]) contains zero-energy modes with 
n = which can be conveniently described in terms of 
Majorana self-conjugated fermionic fieldi^. The ground 
state of the system with multiple spatially well sepa- 
rated vortices with zero bound fermionic modes is topo- 
logically degenerate. The non-Abelyan statistics of vor- 
tex anyons allows for the unitary transformations of the 
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ground state realized through the adiabatic permutation 
of vortices. Such possibility provides an extra motiva- 
tion for the study of vortices in p x + ip y superconductors 
due to their potential application in topological quantum 
computing 1 ^. 

Vortex core Majorana fermions have an important 
property of being stable with respect to the impurity 
scattering^ and order parameter perturbations^. How- 
ever the spectrum vortex core states is extremely sen- 
sitive to the intervortex quasiparticlc tunnelling. The 
corresponding spectrum modification in finite clusters 
of vortices was investigated first by Mel'nikov and 
Silae v 14 ' 15 . In particular for the generic problem of two 
vortices placed at the distance d the low energy fermionic 
spectrum has the form 



£ = — 



± arccos(\/l — e a sin /?) + 7r(n 



1/2 + 7) 
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ip y wave 
= kpd + 
is Gamma 



where 7 = 1/2 for s wave and 7 = for p x - 
superconductivity, a = "lire 
aln(fcF£ 2 /d) + arg(r(l-za)) + 7r/4 where T(x) 
function. The spectrum ^ contains two series of levels 
with the interlevel distance uj ~ Ao/(kp^), where Ao is 
the gap value far from the vortex core, £ is the supercon- 
ducting coherence length. The Eq. @ demonstrates that 
the perturbation of energy levels due to the intervortex 
quasiparticlc tunnelling in general is not described by the 
plain tight binding theory. Indeed the shift of energy lev- 
els with respect to the isolated vortex spectrum becomes 
larger than the energy level spacing u> when the intervor- 
tex distance is smaller than the critical one d < d Cl where 
d c — (£/2)ln(fc^£) is much larger then £ since kp£ ^> 1. 
For the typical parameter kp£ = 100 corresponds to mag- 
netic fields being larger than that of the order 0.lH c2 . 

Thus in a pair of vortices the intervortex quasiparti- 
cle tunnelling removes the twofold degeneracy of vortex 
core states. In particular in p x + ip y superconductor it 
splits the Majorana zero energy states provided the phase 
j3 7^ irn in the Eq.®. Such splitting of zero energy levels 
opens the gap in the fermionic spectrum and can break 
the quantum coherence during the vortex permutation 
which is important for the fault tolerance of topological 
quantum computations^. The generalization of Eq.Q 
for M— vortex clusters is straightforward and was dis- 
cussed in detai l 14 ' 15 . In particular for an odd number of 
vortices M there is always at least one zero energy state 
irrespective of the vortex position in the cluster. Thus 
the splitting of Majorana fermions in vortex cluster is not 
a generic effect and depends on the parity of the total 
number of vortices M. On the other hand in clean type- 
II superconductors without disorder and pinning centers 
vortices in finite magnetic field form periodic Abrikosov 
flux lattice. Therefore the natural question considered in 
the present paper is whether the spectrum of fermions 
on the vortex lattice in chiral p x + ip y superconductor is 
gapped or contains zero energy Majorana states. 

Previously the various types of two dimensional lat- 
tice spectrum problems of Majorana fermions were 



considered ^ 20 ' 22 . These models take into account only 
the tunnelling between lowest energy states in the vortex 
cores. As we have discussed above in the generic problem 
of two vortices the shift of vortex core energy levels be- 
comes larger than the interlevel energy already at small 
magnetic fields H > Q.1H C 2- Therefore the one level 
approximation of lattice models is applicable for sparse 
vortex lattices. Instead in the present paper we consider 
the eigenvalue problem of genuine Bogolubov - dc Gcnncs 
equation in chiral p x + ip y superconductor with gap po- 
tential corresponding to the periodic Abrikosov flux lat- 
tice. To treat this problem we generalize the original ap- 
proach developed earlie r 14 ' 15 to calculate the fermionic 
spectra of finite vortex clusters. This approach allows to 
calculate the spectrum when the intervortex distance is 
d>t 

In general the problem of identifying the quasiparticlc 
energies in superconductors is to solve the Bogolubov - 
de Gennes (BdG) equations having the form: 




(3) 



where H = f 3 [(p - f 3 A) 2 - k 2 F ] /2m, * = (U,V), U 
and V are the particle- and hole - like parts of the 
fermionic quasiparticlc wave function, tr, arc Pauli ma- 
trices, p = — iV, fi are the Pauli matrices in a particle- 
hole space, the gap operator is A = { A(f ), e lx0p } where 
X = ±1 is chirality, f is a coordinate operator, A(r) 
describes the spatial dependence of the gap function and 
{A, B} = (AB + BA)/2 is an anticommutator which pro- 
vides the gauge invariance of A. The phase of the order 
parameter depends on the direction of the electron mo- 
mentum in xy plane: p = p(cos 8 p , sin 9 p ). The magnetic 
field is directed along the z axis B = Bz and for extreme 
type-II superconductors we can consider the magnetic 
field to be homogeneous on the spatial scale of intervor- 
tex distance and take the gauge A = [B x r]/2 where B 
is an average magnetic field. 

Then the periodicity of vortex lattice is determined 
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(4) 
(5) 



where d = n a a + n^b is the translation of the vortex lat- 
tice, n ai rib are integer numbers and ip& is an arbitrary 
constant phase shift. By choosing the Wigner-Seitz ele- 
mentary cell of vortex lattice and placing the origin r = 
at the vortex center in this cell we immediately obtain 
that (p a = iph = tt. 

The translational properties of A(r) and A(r) make 
the Eq.Q to commute with the magnetic translation op- 
erator 



n = T 3 e 



if 3 [Bxd]r/2r 



(G) 



so that T$HT% 



H where Td is the usual translation 
by the lattice vector d. Consequently the solutions of the 
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BdG Eq.Q can be classified according to the eigenstates 
of magnetic translation operator. An important point is 
that the magnetic flux through the vortex lattice unit cell 
is one half of the flux quantum Bz ■ [a x b] = it so that 
the magnetic translations by lattice vectors anticommute 
T a h T£ = -T£t£. Therefore we should introduce the unit 
cell for the quasiparticle functions consisting of two vor- 
tex lattice unit cells, for example shifted by the vector 
a. For the case of square vortex lattice this choice is il- 
lustrated in the FigfTJa). Then the magnetic translation 
subgroup is formed by vectors d m = 2n a a + n&b and the 
solution of Eq. ([3]) in general has the form 



* = J2 e 4Kdm T^ m [l-i (r) + e lKa * 2 (r)] 



(7) 



where the functions ^1,2(1") are localized in the centers 
of vortices forming the unit lattice cell for the quasipar- 
ticles (see FigfT]). We substitute the ansatz (|7|) to the 
BdG Eq.© and calculate the inner product with ^i, 2 (r) 
taking into account only overlap with neighbor cites to 
obtain the system of tight binding equations. 

The further calculation requires expansion of the node 
wave functions ^1,2(1") by the basis of localized fcrmionic 
states of an isolated vortex. It can be implemented us- 
ing the quasiclassical approximation and the convenient 
formalis m 14 ! 15 of the so called s — 8 P representation which 
allows to express the quasiparticle wave function in mo- 
mentum representation in the form: 



+00 



kf J 

—00 

The equation for ip{s,8 p ) reads: Hip = Eip, where 

'o A\ 



(8) 



H = —iVpf 3 d s 



A+ 



(9) 



Vf = kp/m is Fermi velocity. Here we take into ac- 
count the quantization of angular motion variable by 
treating the angular momentum as differential operator 
fi = —id/dOp. Hence the spatial coordinate in the gap 
operator in Eq.Q is quantum variable in s — 6 p repre- 
sentation f = s'kp/kp + {[ki?,z],/i} /k F . Let us empha- 
size that the Hamiltonian © takes account of noncom- 
mutability of jl and 6 P and, thus, the above description 
involves the angular momentum quantization. Replacing 
/t by a classical variable we get Andreev equations along 
straight trajectories. The inner product can be expressed 
through the envelope functions 



/OO />2lT 
ds / d$ P i>tM^p) 
00 Jo 



(10) 



and the magnetic translation operator (j6)) in s — 8 p 
representation has the form T% = f 3 e- iVh ^ 3 / 2 T d 
where <£>h(d) = d(s + nd)sin(6L — 84), the angle 84 de- 
fines the direction of d = d(cos 9^, sin 84) and Xd = 
exp [— idkp(l — ikp l d s )\ is the translation operator—. 



The form of the node functions 'I' 1,2(1") in Eq.([7|) is 
determined by the states localized in isolated vortex. 
We consider the vortex lattice cite 1 centered at the 
origin r = and define the spatial dependence of gap 
function inside the unit cell as A(r) = A v (r)e l9 where 
A„(r = 0) = 0. Then cigenfunctions of the Hamiltonian 
© centered at lattice cites 1 and 2 have the form 



Ms,Op) = C 2 (6 p )T^ v {s,0 P ) 
where n = kp/kp and 

' 1 \ e~ K ^ 



ip v {s,8 p ) 
K(s) = Vf 1 



ir 3 ( 1 +X)8 P /2 



VA 



(11) 
(12) 



(13) 



and A is normalizing factor so 



/ A v (t)dt 


that (vPl^) = 1 ■ The node functions are translated to 
other sides according to the Eq.© so that with the help 
of Eq. (|10[) wc obtain the inner products, e.g. 



(*i|ff|*i>=w 



(* j \H\T£* j ) = i(-iy+ l j b [ 
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(15) 



j = 1,2 where d = a, b and the Hamiltonian H is given 
by ©. The sign in Eq. (|15[) determined by the magnetic 
flux through the unit cell. Here for simplicity we take into 
account the overlap with four neighboring vortices. The 
cases of more neighbors can be considered analogously. 
The main contribution to the inner products of the form 
(|T5|) comes from the stationary points of the phases (akp) 
and (bkp) that is 8* = 8 a b + irn. The stationary points 
8* correspond to the trajectories passing through both 
of the neighbor vortex cores which means that we can 
calculate the overlap factors as follows 



Jd 



7T 

kF J-r 



A(s)-A„(si) ^+{s)f 2 ^ v ( Sl )ds. (16) 



where s\ = s — nd, A„(s) = A v (s)sign(s), and ip v (s) = 
ip v (s,6 p = 0). Then with good accuracy Eq.(fT6f yields 
an estimation |Jd| « A exp(— d/£). 

With the help of the inner products (|14I15[) and taking 
into account that Jd = — J-d we obtain the equations 



(e - wfid = F h (d p )d + F a {8 p )C 2 
(e - ufr)C 2 = -F h (8 p )C 2 + F & {8 P )C X 



(17) 



where F d {8 p ) = Jd sin[d(kF + K)] and d = a, b. The 
system (fT7| should be solved together with periodic 
boundary conditions C\_ 2 {d p ) = Ci^ 2 (8 p + 2tt) tor p x + ip y 
wave and Ci i2 (8 p ) = Ci i2 (6 p + 2ir) for s wave. Note that 
in Eqs. (|17[) we can take into account overlapping with 
next-to neighbor vortices which will introduce the correc- 
tions of the relative order e~ d ^ to the coefficients F a , b- 
Here we neglect such corrections. 
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To solve the Eq. (fl7|) we use the approximate method 
employed earlier for the system of two vortices^. That 
is besides the vicinity of the angles 9* = 9 a: b + irn the 
solution with good accuracy is d 2 - e ^p/". In the 
5 vicinity of the angles 9* the system (fTT)) diagonalizes 
yielding the matching conditions C(9* + S) = MC(6* — 6) 
for the vector C = (C\, C-2) T . The matching matrices are 

M (e b ) = exp(-if 3 X6+) (18) 

M (9 b + tt) = exp(-if 3 X6_) (19) 

M (0 a ) = cosxa+ - in sinxa+ (20) 

M (6 a +7r) = cosx _ -ifisinx a _ (21) 

where Xd± = ^d sin(Kd ± x) where x = ~~ 7r /4 and 
= {Jd/u)^Jir/k F d. 

The periodic boundary condition and Eqs. ([18[) yield 
the Bloch waves e = e„(K) in a periodic Abrikosov flux 
lattice 

/■^ ( arccosX \ . , 

e n ±(K)=u(± ^ +n + 7 j (22) 

where 7 = 1/2 for s wave and 7 = for p x + 
ip y wave, X = cos \a- cos x Q + cos(xb+ + Xb-) - 
sinx a - sinxa+ cos(xt+ - X6-) and n is integer. The 
width of Bloch bands ([22j) is determined by the overall 
amplitude max (to, Aoe~ d ^ / y/kpd) and rapid oscillations 
with the period kp by the intervortex distances a, b. The 
phase of oscillations is determined by the average mag- 
netic field, e.g. a = b = ^/$ /B for the square lattice 
where $0 is magnetic flux quantum for Cooper pairs. 
Notwithstanding the rapid oscillations of the bandwidth 
the spectrum (|22|) survives fluctuations of the vortex po- 
sitions <5r provided their amplitude is relatively small. In- 
deed in case of the disordered vortex lattices let us search 
the quasiparticle waves in the form Vl'(r) = \I/o( r ) + ^( r ) 
where the first term is periodical and given by Eq. ([7]) and 
the second term is the distortion due to the fluctuation 
of vortex positions. We use the expansion (|8I11[) with the 
coefficients Ci,2, (8 P ) to represent the distortion \P(r) at 
the particular lattice site. Then we get the equations for 
Ci,2(# P ) (e — u>ji)Cj = uj£iCjo(9p) where p, = z[Sr x k F ] 
and Cjo(0 p ) corresponds to the periodical part of func- 
tion \l/o ( r ) determined by the solution of Eq. ((TT)) . The 
functions Cjo{9 p ) are the rapidly oscillating ones with 
the characteristic period A9 P rj e d '^/(kF^) so that the 
amplitude of the distortion is small Cj ~ e d ^\5v\/^ <C 1 
provided vortex position fluctuations are small enough 
\St\ < £e~ d/5 - 

The plot of the Bloch band £o+(K) is shown in the 
FigOJb) for the parameters x = 7r /4 and I a = Ih = 0.14. 
One can see that this band contains small energy gap. 
Indeed for |J d | < 1 the Eq.([22]) can be simplified. Tak- 
ing into account the quadratic terms of the order 1^ we 
obtain the gapless spectrum identical to the one level 
lattice modeU^ with the hopping amplitude determined 
by a; I Id |- On the other hand the terms of the order l\ 



open the gap in the spectrum which is beyond the ac- 
curacy of one level approximation and appear due to 
the mixing with higher levels. Note that overlap with 
next-to neighbor vortices introduce correction smaller by 
the factor l/y/kp^ than the interaction with higher lev- 
els. Decreasing the intervortex distance one finally ob- 
tains gapless spectrum when |id| > tt. The example 
of such crossover to the gapless regime is demonstrated 
in the dependence e(K x ,K y ) for the square lattice in 
FiglHc) . We set x = 7r /4 and plot by red dash-dotted, 
blue dashed and black solid lines the Bloch bands £o± (K) 
for J a = Z a = 1.4; 2.8; \/2ir correspondingly. For higher 
values of the overlap |/d| the structure of Bloch bands be- 
comes more complicated with rapid oscillations as func- 
tion of quasimomentum with the characteristic period of 
the order l/(\Id\d). Such complicated structure of Bloch 
band £o+(K) for x = t/4 and I a = I a = 4.3 is shown in 
the contour plot FigJIJd). 

Finally let us consider the possible experimental 
test of the suggested gapless spectrum of Majorana 
fermions (f!?2")) . 

In addition to the variety of experiments proposed^, 
the electronic thermal conductivity n measurements have 
been proven as an effective tool to study the quaiparti- 
cle spectrum in the vortex phased. The electronic states 
in magnetic Bloch bands (|22|) can carry the energy cur- 
rent in the direction _L B due to the hopping of quasi- 
particles between neighboring vortices. The electronic 
spectrum on Abrikosov lattice (|22|) is gapless even in 
the regular lattices provided the intervortex distance is 
d < d c . Hence one should expect the threshold behav- 
ior of n{B) in the increasing magnetic field in the limit 
T — > 0. That is k should be zero in the gapped phase 
and in the gapless regime it can be estimated by the 
textbook expression n± ~ TVgTv where V g = de/dK. 
is the group velocity of Bloch waves (f2"2"j) . r is transport 
time and v = cu~ 1 d~ 2 is density of in the vortex lattice. 
The typical value of group velocity determined by the 
intervortex hopping is V g ~ Vpe~ d ^ yjd/kp^ 2 so that 
assuming d = /3£ \J H C 2 / B where /3 ~ 1 we obtain 



^/^ = ^^e-^V^. (23) 

Besides that n± also contains oscillating part due to rapid 
oscillations of energy levels (|22[) with the period kp 1 by 
the intervortex distance. However the oscillations are 
mostly cancelled out due to the complicated structure of 
Bloch bands shown in FigQJd). The obtained value of k 
is valid in the superclean limit lot > 1 for fully gapped 
superconductors including s and p x + ip y wave symme- 
tries. The estimation ([23| contains small prefactor 
which can explain the experimentally observed small val- 
ues of k± at B <C H c -2r^. Interestingly the observed^ 
similar to (|23|) behavior of K|| in the direction || B can be 
explained by the theory-^ applied to the spectrum ((2"2"j) 
since the number of conducting modes along vortex line 
is determined by the tunnelling factors Id- 
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